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The  term  “water  boiler"  reactor  refers  to  a  type  of  aqueous  homogeneous  reactor  (AHR)  that  was 
designed,  built  and  operated  by  Los  Alamos  in  the  1940s.  This  was  the  first  type  of  liquid  fuelled  reactor 
and  the  first  to  be  fuelled  with  enriched  Uranium.  For  security  reasons  the  term  “water  boiler"  was 
adopted  and  three  versions  were  built:  LOPO  (for  low  power),  HYPO  (for  high  power)  and  SUPO  (for  super 
power)  which  were  spherical  shaped  reactor  vessels.  The  name  was  appropriate  as  the  reactors  appeared 
to  boil  although  this  was  actually  due  to  the  release  of  radiolytic  gas  bubbles;  although  SUPO  was  oper¬ 
ated  during  some  studies  close  to  the  boiling  point  of  uranyl  nitrate.  The  final  water  boiler  “SUPO"  was 
operated  almost  daily  as  a  neutron  source  from  1951  until  its  deactivation  in  1974-23  years  of  safe,  reli¬ 
able  operation.  Many  of  the  key  neutron  measurements  needed  in  the  design  of  the  early  atomic  weapons 
were  made  using  LOPO,  HYPO  and  SUPO.  More  recently  SUPO  has  been  considered  as  a  benchmark  for 
quasi-steady-state  operation  of  AHRs  with  internal  cooling  structures. 

This  paper  presents  modelling  and  analysis  of  the  coupled  neutronic  and  fluid  time  dependent  charac¬ 
teristics  of  the  SUPO  reactor.  In  particular  the  quasi-steady-state  dynamics  of  SUPO  have  been  investi¬ 
gated  together  with  its  heat  transfer  characteristics.  In  the  simulations  presented  the  SUPO  reactor  is 
modelled  using  the  spatially  dependent  neutron/multiphase  CFD  simulation  tool,  FETCH,  at  a  quasi¬ 
steady-state  power  of  25  kW.  SUPO  also  possessed  a  cooling  coil  system  that  fed  cooling  water  through 
the  reactor  for  the  extraction  of  the  fission  and  decay  heat.  This  cooling  system,  and  the  heat  extraction,  is 
modelled  in  the  simulations  using  a  new  sub-modelling  approach  that  is  detailed  here.  The  results  from 
this  simulation,  such  as  gas  fraction,  gas  generation  rate,  coolant  rate  and  average  temperature,  are  com¬ 
pared  against  the  available  experimental  information. 

©  2012  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

The  development  of  the  water  boiler  reactor  at  Los  Alamos 
(Bunker,  1963;  King,  1956, 1990;  Los  Alamos  Scientific  Laboratory, 
1951;  Bunker,  1983;  Rosenthal,  2003;  King,  1952;  Beck  et  al„  1952; 
Lyon,  1953;  Thomas,  1953;  Kasten,  1954;  Shorthall  et  al„  1954; 
North  American  Aviation,  1954;  Durham,  1955;  Nelson  and  Mann, 
1955;  Lane  et  al„  1958;  Gamble,  1959;  Stitt,  1959;  Flora  et  al„ 
1962;  Dunenfeld,  1962;  Spiegler  et  al„  1962;  Schulberg,  1965; 
IAEA,  2008),  was  made  possible  by  the  first  availability  of  enriched 
uranium  from  Oak  Ridge,  and  this  acted  as  a  supporting  demon¬ 
stration  to  the  weapons  programme  on  building  and  modelling 
critical  assemblies.  Modelling  had  suggested  that,  with  a  suitable 
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reflector,  criticality  could  be  achieved  with  a  fuel  in  the  form  of  a 
uranyl  sulphate  solution  when  contained  within  a  spherical  stain¬ 
less  steel  tank  of  diameter  approximately  30  cm.  It  was  decided  to 
construct  three  reactors  over  5  years.  The  first  was  a  low  power 
uranyl  sulphate  LOPO  reactor  and  this  was  followed  by  the  high 
powered  (HIPO)  and  super  powered  (SUPO)  reactors  which  con¬ 
tained  a  uranyl  nitrate  fuel.  The  SUPO  design  was  developed  in 
1951  (Bunker,  1983;  Rosenthal,  2003)  and  operated  with  a  fuel 
that  was  enriched  to  88%  and  was  surrounded  by  an  all  graphite 
reflector.  The  highly  enriched  fuel  meant  that  the  reactor  could 
reduce  the  amount  of  nitric  acid  decomposition  whilst  enhancing 
its  reactivity.  The  reactor  was  safely  operated  for  23  years  before 
being  decommissioned,  and  was  regularly  run  at  a  powers  of 
25  kW  or  more  for  durations  that  lasted  several  hours  at  a  time. 
The  success  of  these  early  water  boilers  resulted  in  the  reactor 
technology  spreading  throughout  the  US  and  around  the  world  to 
countries  including  Denmark,  Japan,  Germany  and  Italy. 
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Pictures  of  the  SUPO  reactor  taken  at  the  time  of  its  construction 
are  presented  in  Figs.  1  and  2.  These  show  the  reactor  core  vessel  to 
have  a  spherical  geometry  and  the  coils  within  its  body  are  for  the 
feeding  of  cooling  water  for  the  extraction  of  fission  and  decay 
heat.  Three  separate  sets  of  cooling  coils  were  used  which  allowed 
higher  operating  powers  to  be  studied.  The  SUPO’s  main  reactivity 
control  was  achieved  through  two  boron  control  rods  moving 
within  vertically  aligned  guide  tubes  -  see  Fig.  1  (bottom).  The  ves¬ 
sel  also  possessed  a  ‘glory  hole’  (clearly  shown  in  Fig.  2)  that 
passed  through  its  near  central  region  and  which  was  used  for 
the  production  of  neutron  beams. 

The  SUPO  reactor’s  design  and  composition  has  a  number  of 
interesting  aspects  when  one  considers  the  various  physical  pro¬ 
cesses  that  determine  its  behaviour  and  ability  to  operate  over 
periods  lasting  several  hours.  Most  notable  would  be  the  combina¬ 
tion  of  the  neutronic  and  fluid  properties  together  with  the  various 
influences  they  have  on  each  other.  In  particular,  radiolytic  gases 
that  are  produced  by  fission  affect  both  the  neutronic  population 
and  fluid  dynamic  flow  patterns.  That  is,  bubbles  that  form  within 
the  fuel  solution  increase  voidage,  thus  lowering  reactivity,  whilst 
their  migration  to  the  solution  surface  impose  drag  forces  that 
influence  the  fluid  flow.  Other  properties  that  couple  the  neu¬ 
tron-fluid  physics  include  the  migration  of  the  delayed  neutron 
precursors  that  flow  within  the  fluid,  fission  heat  sources  that  af¬ 
fect  fluid  density  and  flow,  and  variations  in  the  nuclear  cross-sec¬ 
tions  with  respect  to  fluid  density  and  temperature.  These 
processes  are  further  affected  by  the  SUPO  geometry  and  internal 
cooling  coil  and  guide  tube  structures.  The  presence  of  these  will 
clearly  affect  the  flow  of  the  fluid  due  to  drag,  but  also  the  heat 


Fig.  2.  View  of  vessel’s  internal  parts  including  the  inner  (traced  in  red)  and  outer 
(traced  in  blue)  cooling  coils,  glory  hole  and  control  rod  guide  tubes  (Bunker,  1963; 
King,  1 990).  (For  interpretation  of  the  references  to  colour  in  this  figure  legend,  the 
reader  is  referred  to  the  web  version  of  this  article.) 

exchanges  with  the  cooling  coils  will  further  affect  the  solution 
(and  neutron)  dynamics.  The  SUPO  reactor  therefore  provides  a  un¬ 
ique  benchmark  for  modelling  solution  reactors  and  other  AHRs, 
particularly  with  imbedded  cooling  coil  systems.  It  is  therefore  also 


Fig.  1.  Cooling  coils  viewed  from  below  the  upper  vessel  hemisphere.  Above  left;  all  three  cooling  tubes  in  place  plus  location  of  CR  guide  tubes  and  glory  hole  (Bunker,  1 963, 
King,  1990).  Above  right:  the  outer  cooling  coil  (King,  1990).  Bottom:  all  tubes  seen  from  axis  with  CR  guide  tubes  present  (King,  1990). 
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an  ideal  system  for  model  verification  for  determining  the  various 
physical  properties  of  reactors  containing  uranyl  nitrate  fuels.  Key 
parameters  such  as  the  rates  of  radiolytic  gas  production  and  heat 
production/extraction  as  well  as  knowledge  of  the  fissile  fluid 
flows  all  play  vital  roles  in  the  safe  and  efficient  operation  of  such 
reactor  types.  The  importance  to  be  able  to  accurately  model  such 
multi-physics  phenomena  is  therefore  clear  to  be  seen. 

This  article  presents  a  study  on  the  time-dependent  and  quasi- 
steady-state  behaviour  of  the  neutron-fluid  dynamics  of  the  SUPO 
reactor  by  its  simulation  using  the  FETCH  model  Pain  et  al.  (2001a). 
FETCH  is  a  multi-physics  code  that  solves  time-dependent,  coupled 
neutron  transport  and  fluid  dynamic  problems  for  nuclear  systems, 
and  this  includes  those  involving  fissile  solutions.  The  temporal 
and  spatial  dependence  of  the  neutron  field  is  resolved  using  the 
EVENT  module,  de  Oliveira  (1986),  which  is  based  on  the  second 
order  even-parity  formulation  of  the  Boltzmann  transport  equa¬ 
tion.  The  fluid  dynamics  is  resolved  in  the  context  of  multi-phase 
phenomena  and  this  is  through  the  module  FLUIDITY  Pain  et  al. 
(2001b).  FETCH  provides  the  interface  between  these  two  modules. 
It  couples  them  by  resolving  the  various  dependencies  the  two 
physical  fields  have  on  each  other,  and  this  is  shown  in  Fig.  3.  This 
illustrates  the  heat  sources,  fission  distribution  and  delayed  neu¬ 
tron  precursors,  that  are  created  within  the  neutronics  component 
through  fission,  being  fed  into  the  fluid  calculation.  In  return,  the 
spatial  distribution  of  the  delayed  neutron  precursor,  material 
temperature,  density  and  void  fractions  (radiolytic  gas  to  liquid 
volume  ratios)  are  fed  back  to  the  neutronics  calculation.  These 
feedback  effects  enable  the  detailed  modelling  of  fissile  solutions 
such  that,  for  example,  the  production,  transport  and  decay  of 
the  delayed  neutron  precursors  can  be  tracked  whilst  they  flow 
within  the  fluid.  In  addition,  fission  rates  are  used  to  calculate  both 
dissolved  radiolytic  gas  production  and  bubble  formation.  These 
bubbles  are  represented  within  the  multi-phase  context  through 
a  gas  phase  that  is  modelled  together  with  the  liquid  phase  of 
the  fuel.  This  enables  the  bubbles'  migration  towards  to  the  solu¬ 
tion’s  surface  to  be  included  in  the  simulation  and  so  their  appear¬ 
ance  and  disappearance  will  be  reflected  in  the  power  profile.  The 
full  coupling  of  these  components  in  FETCH  has  helped  understand 
the  dynamics  of  the  SUPO  reactor  and  provide  power  traces  that 
fluctuate  subject  to  all  these  variables.  These  aspects  can  only  be 


fully  understood,  and  therefore  modelled,  by  considering  both 
the  neutronics  and  fluid  properties. 

In  the  FETCH  model,  both  the  EVENT  and  FLUIDITY  components 
apply  finite  elements  in  resolving  the  spatial  dependence  of  their 
respective  fields.  They  can  resolve  geometries  represented  in  both 
3D  Cartesian  coordinates  and  2D  axisymmetric  (RZ)  domains.  The 
details  of  the  two  models  are  briefly  reviewed  here,  and  this  in¬ 
cludes  a  description  of  the  equations  describing  radiolytic  gas  bub¬ 
ble  production,  transport  of  delayed  neutron  precursors,  fission 
heat  transfer  onto  the  fluids  and  macroscopic  cross-section  gener¬ 
ation  (which  vary  with  respect  to  changes  in  fluid  density,  due  to 
gas  production,  and  temperature).  The  FETCH  model  also  allows 
for  the  movement  of  the  control  rod  during  a  simulation  and  in¬ 
cludes  an  algorithm  that  adjusts  its  height  in  order  to  achieve 
and  maintain  a  desired  power.  In  addition  to  this,  in  order  to  sim¬ 
ulate  the  heat  loss  from  the  fissile  solution  to  the  cooling  coil  sys¬ 
tem,  a  new  sub-model  has  been  developed.  This  uses  separate  1 D 
models  to  represent  each  of  the  cooling  coil  pipes  and  these  ex¬ 
change  heat  with  the  main  model’s  fissile  solution  and  advects  this 
through  the  pipe  system.  The  processes  of  heat  exchanges  between 
the  fissile  solution  and  pipe  model  are  described  in  full  and  com¬ 
parisons  of  the  modelled  outlet  temperature  are  made  against 
SUPO  data.  The  simulation  of  the  SUPO  reactor  is  based  on  it  in 
operation  at  25  kW  of  power  which  will  achieve  a  temperature  of 
approximately  60  °C.  The  control  rod  positions  are  adjusted  in 
the  model  so  that  they  respond  to  changes  in  the  system’s  reactiv¬ 
ity  and  so  maintain  this  intended  operating  power.  The  influence  of 
the  speed  of  the  control  rod’s  movement  on  the  reactor’s  behaviour 
is  studied.  This  is  accompanied  by  the  predicted  power  levels,  solu¬ 
tion  pressures  and  velocities,  as  well  as  the  outflow  temperatures 
of  the  cooling  coil  water.  The  results  from  these  simulations  are 
compared  against  the  available  experimental  information  (Bunker, 
1963,  King,  1956,  1990). 

The  following  sections  are  set  out  as  follows.  Section  2  describes 
the  SUPO  design  and  fuel  specifications  with  the  reactor  operating 
at  its  intended  25kW  power  level.  This  includes  the  description  of 
the  FETCH  model  that  recasts  the  reactor  geometry  in  RZ  dimen¬ 
sions.  The  radiation  and  multi-phase  flow  equations  solved  by 
the  FETCH  model  are  then  described  in  Section  3,  and  this  includes 
a  description  of  the  coupling  between  the  physical  fields.  The 
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description  of  the  cooling  coil  sub-model,  that  models  the  heat 
extraction  system,  has  been  presented  in  its  own  section  -  Section 
4.  Section  5  presents  the  results  of  the  FETCH  simulations  and  the 
conclusion  are  presented  in  Section  6. 

2.  The  SUPO  design  and  specifications 

The  SUPO  reactor  design  and  means  of  operation  drew  heavily 
on  experience  with  LOPO  and  HYPO.  The  reactor  vessel  described 
in  Section  1  was  enclosed  and  centred  within  a  graphite  reflector 
which  was  in  turn  extended  in  the  form  of  two  graphite  thermal 
columns.  This  is  illustrated  in  Fig.  4  which  shows  a  cut-away  draw¬ 
ing  of  the  complete  SUPO  reactor.  This  graphite  was  represented,  in 
simplified  geometry,  in  the  EVENT  neutronics  component  of  the 
coupled  FETCH  calculations.  Concrete  shielding  enclosed  the 
reflector  and  thermal  columns.  In  addition  to  solution  height  mea¬ 
suring  devices  within  the  vessel  there  were  recombiner  provisions 
(seen  mounted  above  the  vessel  in  Fig.  4)  to  capture  fission  gases 
and  gaseous  products  of  radiolytic  decomposition  and  also  to  re¬ 
turn  to  the  vessel  aqueous  condensate.  The  provision  of  guide 
tubes  for  two  in-vessel  boron  control  rods  has  been  referred  to; 
the  boron  was  enriched  in  boron-10.  In  addition,  two  cadmium 
metal  strip  rods,  very  likely  used  as  safety  rods,  moved  in  vertical 
aluminium  frames  within  the  graphite  and  tangential  to  the  vessel 
outer  surface.  The  drives  for  the  boron  rods  and  the  aluminium 
frames  for  the  cadmium  rods  may  be  seen  in  Fig.  4.  In  her  report, 
Bunker  (1963)  describes  the  25  kW  reference  case  and  also  gives 
control  rod  worth  data.  In  addition  to  the  glory  hole,  a  second  tube 
tangential  to  the  vessel  exterior  may  be  seen.  In  order  to  reduce 
peak  solution  temperatures  within  the  vessel,  presumably  to  re¬ 
duce  the  margin  before  boiling,  the  inlet  water  to  all  three  cooling 
tubes  was  reduced  to  5  °C  using  a  chiller.  Bunker  also  discusses  the 
detailed  heat  balance  in  SUPO  under  operating  conditions:  as  a 
proportion  of  total  fission  power  she  ascribes  ~2%  and  ~6%, 
respectively,  to  heat  removed  by  the  condenser  and  to  heat  lost 
by  the  vessel  through  conduction  into  the  graphite  reflector.  The 
present  authors  expect  to  perform  modelling  to  investigate  these 
factors  in  future  work.  The  factors  given  by  Bunker  are  broadly 
consistent  with  data  given  by  Durham  (1955). 


Fig.  4.  Cutaway  drawing  of  SUPO  (Bunker,  1963;  King,  1990) 


As  described  by  Bunker,  the  SUPO  reactor  has  a  spherical  shell 
made  from  two  welded  hemispheres  constructed  from  stainless 
steel  with  thickness  0.15  cm  and  outer  dimension  30.48  cm.  It  con¬ 
tains  a  soup  of  uranyl  nitrate  that  occupies  a  volume  of  12,572  cm3 
(at  60  °C)  with  Table  1  listing  a  breakdown  of  its  constituent  iso¬ 
topes.  The  three  cooling  coils  have  identical  submerged  length  of 
20  ft,  they  have  inner  and  outer  diameters  0.476  cm  and 
0.635  cm  respectively,  and  displace  the  fluids  by  a  total  volume 
of  550  cm3.  The  control  rod  guide  thimbles  enter  the  solution  to 
a  depth  of  approximately  23  cm,  with  each  displacing  the  fluids 
by  64.9  cm3.  The  glory  hole  is  a  stainless  steel  tube  of  diameter 
2.8  cm,  and  enters  the  vessel  by  a  length  of  27.28  cm  which  dis¬ 
places  the  fluids  by  240  cm3.  These  specifications  are  considered 
the  basic  data  set  describing  the  conditions  of  the  reactor. 


2.1.  The  RZ  semi-explicit  model  of  the  SUPO  reactor 

The  FETCH  model  used  in  the  SUPO  simulation  has  made  use  of 
the  reactor’s  spherical  geometric  features  and  has  been  repre¬ 
sented  in  2D  RZ  geometry.  This  was  chosen  as  it  allows  the  use 
of  a  fine  spatial  resolution  mesh  that  can  resolve  much  of  the 
dynamics  around  the  cooling  coils  whilst  keeping  CPU  run  time  rel¬ 
atively  low.  The  cooling  coils  within  the  FETCH  model  are  explicitly 
resolved,  and  as  a  result  much  of  the  physics  behind  SUPO  flow 
dynamics  can  be  determined.  Although  the  3D  dynamics  of  SUPO 
may  be  different  in  detail  from  the  RZ  dynamics,  the  2D  approach 
can  (in  the  absence  of  a  parallel  highly  resolved  3D  FETCH  model) 
yield  results  distinct  from  those  from  a  3D  model.  This  is  because  a 
3D  simulation  will  require  a  more  parameterised  model  when  run¬ 
ning  in  a  serial  environment. 

The  key  geometric  features  of  SUPO  requires  an  unstructured  fi¬ 
nite  element  mesh  to  accurately  represent  the  solution  geometry 
as  well  as  the  surrounding  steel  vessel  and  graphite  (see  Fig.  5). 
Each  intersection  of  the  cooling  coil  with  the  RZ-geometry  is 
explicitly  resolved.  The  volume  of  each  intersection  is  maintained 
so  that  it  is  the  same  as  the  cooling  coil  pipes  cross  sections.  To 
simplify  the  explicit  meshing  of  cooling  coils  pipes,  their  geometry 
was  made  rectangular  with  the  same  equivalent  cross  sectional 
area.  The  model  mimics  the  idealised  geometry  of  the  pre-con¬ 
struction  engineering  drawings,  and  the  distribution  of  cooling 
coils  are  represented  as  sets  of  nesting  helices.  The  three  coils 
are  represented  by  six  separate  sets  of  vertically  aligned  circular 
rings  of  tubing.  Each  tube  is  609.6  cm  (20  ft)  in  length  and  makes 
up  an  appropriate  pair  (the  cooling  water  flows  down  one  coil 
and  up  another),  together  with  the  necessary  length  of  pipe  to  con¬ 
nect  with  the  vessel  head.  The  dimensions  of  the  cooling  coils  are 
listed  in  Table  2.  Due  to  the  dominance  of  the  circumferential  com¬ 
ponent  of  the  tube  length  for  a  helix,  only  for  the  first  inner  set  of 
tube  rings  is  it  necessary  to  correct  for  the  non-horizontal  align¬ 
ment  of  the  original  helical  structures.  That  is,  the  helices  are 
slightly  adjusted  to  ensure  the  correct  quantity  of  mass  is  sub¬ 
merged  within  the  solution.  Finally,  although  the  glory  hole  is 
completely  ignored  in  this  geometry,  the  control  rod  has  been  rep¬ 
resented  but  its  presence  is  only  made  to  be  felt  by  the  neutrons. 
The  control  rod  guide  tube  can  be  seen  to  enter  the  reactor’s  vessel 
through  its  top  right  section,  and  the  control  rod  can  move  verti¬ 
cally  within  this  sleeve  to  adjust  reactivity.  However  the  sleeve 
itself  is  arranged  so  as  not  to  affect  the  solution  flow  patterns.  This 


volume  of  12,700  cm3. 


Isotope  U235  U238 

Weight  (g)  950  120 


that  forms 


N 

190 


a  solution  with  total  a 


H  O 

1410  11,960 
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a  b 


Fig.  5.  (a)  Geometry  of  SUPO’s  fuel  section  including  the  mesh  and  the  intersection  of  the  cooling  coils  with  the  RZ  domain.  There  are  849  quadrilateral  finite  elements  and 
899  nodes  and  this  is  the  domain  used  for  the  fluids  part  of  the  FETCH  calculation.  A  finer  mesh  is  produced  by  splitting  each  element  into  four  and  thus  there  are  4  x  1267 
elements  in  the  fine  mesh  simulation,  (b)  The  sequence  and  direction  of  flow  of  the  cooling  water  within  the  three  sets  of  cooling  coils  are  outlined  by  the  black  lines.  The 
region  occupied  by  the  control  rod  is  also  indicated  -  although  its  presence  will  only  be  felt  by  the  neutronic  field. 


Table  2 

Dimensions  and  flow  directions  of  the  cooling  coil  helices  used  in  SUPO  model. 


Ring  (Pipe)  and 


1  (l)down 

2  (2)  down 

3  (2)  up 

4  (l)up 

5  (3)  down 

6  (3)  up 


Radius  Circumference  Separate  Tube 


(cm)  (cm) 


1.95  12.25 

4.40  27.65 

6.94  43.60 

8.65  54.35 

11.12  69.87 

13.35  83.88 


rings  length 

(on) 

9  110.3 

9  248.8 

8  348.9 

9  489.4 

5  349.4 

3  251.6 


is  because  its  representation  and  location  in  the  RZ  geometry  will 
create  a  solid  barrier  between  the  inner  and  outer  regions  of  the 
reactor’s  upper  hemisphere.  Instead,  the  fluid  can  flow  through 
the  control  rod  section  as  though  it  was  not  there. 

The  mesh  in  the  SUPO  reactor  has  been  designed  not  to  contain 
a  large  number  of  finite  elements  so  that  the  FETCH  simulations 
can  run  rapidly.  A  FETCH  simulation  of  SUPO  takes  about  2  h  to 
perform  approximately  30  min  of  real  time.  It  takes  about  10  min 
of  simulation  time  for  the  SUPO  system  to  reach  a  quasi-steady- 
state  where,  in  a  time-averaged  sense,  the  heat  losses  (to  the  cool¬ 
ing  coil  water)  are  balanced  by  the  fission  heat  source  and  the  con¬ 
trol  rod  algorithm  is  maintaining  a  quasi-steady-state  condition. 


3.  The  FETCH  model  for  SUPO 

This  section  provides  a  brief  outline  of  the  FETCH  model  and 
states  the  equations  that  are  solved  for  modelling  the  transport 
of  neutrons  and  the  solution’s  fluid  dynamics.  In  modelling  the 
SUPO  reactor  the  solution  is  resolved  in  the  context  of  multi-phase 
flow  where  it  is  considered  as  a  mixture  of  both  liquid  and  gas 
phases.  The  latter  of  these  phases  is  used  to  account  for  the  produc¬ 
tion  of  radiolytic  gas  bubbles  within  the  simulation.  This  section 
also  includes  a  description  of  the  coupling  between  the  radiation 
and  fluid  fields  and  describes  how  fluid  heat  sources  are  generated, 
how  dissolved  gas  and  bubbles  are  produced  and  how  these  feed 
back  into  the  nuclear  material  cross-section  data. 


3.1.  Neutron  transport  module:  EVENT 

EVENT  is  a  general  purpose  radiation  transport  model  used  in 
radiation  shielding,  reactor  and  criticality  calculations  in  both  tran¬ 
sient  and  non  transient  modelling.  Its  formulation  is  based  on  a 
second  order  even-parity  principal,  which  can  be  formed  from 
the  first  order  Boltzmann  transport  equation  that  describes  the 
transport  of  neutral  particle  radiation, 

i  ^  iA(r,  ft,  E,  t)  +  ft  •  ViA(r,  ft,  E,  t)  +  7#(r,  ft,  E,  t) 

« <S(r,  ft,  E,t).  (1) 

The  variable  xj/(r,  £1,  E,  t)  denotes  the  particle  intensity  (or  angu¬ 
lar  flux)  in  the  7  dimensional  phase-space  of  space  (r),  angle  (£1), 
Energy  (E)  and  time  (t).  The  particle  velocity  is  denoted  by  v.  The 
first  term  of  the  equation  describes  the  rate  of  change  in  the  parti¬ 
cle  distribution  (w.r.t.  time)  through  a  small  volume  in  phase- 
space.  This  accounts  for  particles  streaming  out  of  the  volume 
(the  second  term)  as  well  as  those  being  absorbed  and  scattered 
by  the  material  (the  third  term).  This  scattering  and  removal  oper¬ 
ator  is  expressed  fully  as: 


H<j/( r,  ft,  E,  t )  =  fft(r,  E,  ft,  E,  t)  -  ^  dE'  jf  dft'trs(r;  ft' 

— >  ft,£'  — >  E,t)f(r,£i,E!,t).  (2) 

where  nt  and  os  denote  the  macroscopic  cross-sections.  The  term  on 
the  right  of  Eq.  (1)  accounts  for  any  source  of  particles  which  may 
originate  from  external  sources  (5ext),  through  fission  events  or 
from  the  decay  of  delayed  precursors: 

S(r,  ft,  E,  t)  =  <Sext(r,  ft,  E,  t )  +  ^^Zp(£)  j  dft' 

x  £  dE' voy(r,  E',  t)*(r,  ft',  E't, )  +  $>z2(JE)  Dk.  (3) 

The  term  Dk(r,t)  denotes  the  precursor  concentration  for  de¬ 
layed  group  k  and  voy  is  the  average  number  of  neutrons  emitted 
per  fission,  multiplied  by  the  fission  cross-section,  p  is  the  sum  of 
delayed  precursor  fractions  which,  in  this  article,  is  comprised 
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from  6  precursor  groups  that  have  half  lives  Xk.  The  terms  /p(£)  and 
Xi ‘(E)  denote  the  spectrum  of  the  prompt  fission  neutrons  and  de¬ 
layed  fission  neutrons  (from  precursor  group  k)  respectively. 

By  applying  a  time  discretisation  such  backward  Euler  to  the 
time  term  of  Eq.  (1),  one  can  fold  the  temporal  information  within 
the  scatter/removal  and  source  terms  to  form  the  following  expres¬ 
sion  of  the  flux  at  the  end  of  a  time  step: 
fl  •  V iA(r,  n,  E,  t)  +  S1,E,  t)  =  <S*  (r,  SI,  E,  t) .  (4) 

The  terms  introduced  here  are  defined  as  H*  =  H  +  ^  and 
5*  =  S  +  ^  with  i j/0  denoting  the  flux  at  the  beginning  of  the  time 
step.  This  enables  one  to  solve  a  sequence  of  time  independent 
equations  where  the  flux  at  each  time  step  is  found  using  its  solu¬ 
tion  from  the  previous  step  as  a  component  of  its  source.  Now 
using  a  multi-group  approach  for  treating  the  energy  dependence, 
a  mono  energetic  time  independent  transport  equation  can  be  con¬ 
sidered  (i.e.  E  is  now  dropped).  The  solution  procedure  is  to  then 
form  the  even  and  odd  parity  variables  of  the  angular  flux  (t //+ 
and  «/r  respectively): 

r  =  Wr(O)  +  *(-n)]/2  and  r  =  Wr(ll)  -  (5) 

together  with  similar  parity  variables  of  the  source  term, 

s*-.[s(n)±a(-o)]/2.  (6) 

Eq.  (4)  is  then  expressed  in  terms  of  —SI: 

-SI  •  ViKr,  n,  t)  +  H*f( r,  -si,  t)  m  <S*(r,  -si,  t),  (7) 

which  is  then  added  and  subtracted  from  Eq.  (4).  By  substituting  in 
the  expressions  (5)  and  (6),  the  following  two  sets  of  equations  are 
formed: 

Ci/r+=S+-Sl-  Vi/T,  (8) 

(9) 

The  two  operators  C  and  C  define  the  removal  terms  minus  the 
even  and  odd  components  of  the  distributed  scattering  H  operator, 
respectively. 

Ci/f(r, «)  =  fft i/r(r, n)-  J4  dCl'(T+( r,  Sl.Sl')i/r(r,  a'),  (10) 

G“V(r,n)  =  (7t^(r,  SI)  -  J  dSl'a-(r-,Sl.Sl')i//(r,Sl').  (11) 

Eqs.  (8)  and  (9)  can  now  be  combined  to  eliminate  the  odd  parity 
moments  and  form  the  second  order  differential  equation  only  in 
terms  of  the  even  parity  variables. 

-si  ■  vgsi  •  ViA+  +  Cxi/+  =  s+-  sivgs~  (i  2) 

In  obtaining  the  solution  to  this  formulation,  the  neutron  scalar 
flux  at  position  r  can  be  recovered  through  the  following  relation¬ 
ship  with  the  even-parity  function: 

<i !>(r)  =  J  i(r,Sl)dSl  =  J  xlr+(r,Sl)dSl,  (13) 

from  which  the  power  of  the  reactor  at  a  given  time  instance  is  gi¬ 
ven  as: 

Power  =  WJ  J  0/</>(r,E).  (14) 

where  w  represents  the  energy  released  per  fission  -  w  =  200  MeV. 
A  similar  relationship  can  be  defined  relating  the  solution’s  current 
with  the  odd  parity  variables. 

J(r)=  J  Shls(r,Sl)dSl=  J  Shl/~(T,Sl)dSl.  (15) 

The  solution  of  Eq.  (12)  can  be  shown  to  be  equivalent  to  mni- 
mising  the  functional: 


=  (nViA,S«Vi/0  +  -  2(^,5+)  -  2(ftViA,eS“),  (16) 

which  is  the  approach  taken  in  the  EVENT  module,  with  full  details 
given  in  Pain  et  al.  (2001a).  The  discretisation  methods  used  to  re¬ 
solve  the  dimensions  of  this  equation  are  as  follows.  Spherical  har¬ 
monics  are  used  to  discretise  the  angular  dimensions,  and  in  this 
article  a  P,  resolution  is  used  throughout.  As  already  stated,  bi-lin- 
ear  finite  elements  are  used  to  discretise  the  spatial  domain,  Back¬ 
ward  Euler  is  applied  in  time  and  multi-group  in  energy. 

3.2.  Multi-phase  solver:  FLUIDITY 

FLUIDITY  solves  the  multi-phase  flow  equations  which,  in  the 
context  of  this  article,  involves  the  two  phases  of  liquid  and  gas. 
For  each  phase  three  conservation  equations  are  solved  that  con¬ 
serve  the  mass,  momentum  and  thermal  energy.  The  mass  equa¬ 
tions  are  used  to  track  the  volume  fraction  ak,  density  pk  and 
velocity  vk  of  each  phase  k  and  is  given  by  the  equation, 

kPk )  +  V  ■  ( ctkpkvk )  -  V  •  pkyckW ctk  -S^  =  0,  (17) 

where  S™  denotes  a  source  of  mass.  In  the  momentum  equations  the 
momentum  exchanges  between  the  phases  are  conserved.  How¬ 
ever,  they  also  account  for  interfacial  pressure  Mpk  and  drag  forces 
Mk,  virtual  mass  Mkm  and  lift  Mk,  Soo  (1990).  For  phase  k  the 
momentum  equation  is  given  as, 

0y.kpkvk  +  v  akf)kVkVk  _  _akvp  +  V  ■  ak(rk)  +  a kpkg 

+  (Gkuk  -  G^tv)  +Mk  +  Mk 
+  Mkm  +  Mk,  (18) 

where  the  force  terms  are  defined  as  force  per  unit  volume,  p  is  the 
shared  pressure  of  all  phases,  g  the  gravitational  acceleration  and  xk 
the  viscous  stress  tensor.  The  viscous  stress  tensor  is  given  by 
xk  —  pk  1  [ V vk  +  V vk,  where  the  superscript  T  denotes  the  trans¬ 
pose  of  V  vk  and  pk  is  the  viscosity  of  phase  k.  The  terms  involving 
Gk,Gk’  represent  the  momentum  exchange  between  the  phases 
due  to  a  transfer  of  mass  (momentum)  from  one  phase  to  another 
and  visa  versa.  The  subscript  k!  is  attached  to  variables  associated 
with  phases  other  than  those  of  the  liquid  and  gas. 

The  equation  for  conserving  internal  energy  of  the  liquid  phase 
e(  is  given  by: 

J^(a,p,e,)  +  V  ■  (a,p,z/,e,)  +p(V  •  ct,v,  +  V  •  otgvg) 

=  V  ■  a  ,yfVT,  +  Q,  +  r,.  (19) 

It  is  assumed  that  ei  is  related  to  liquid  temperature  Tj  by 
e;  =  CviTi,  where  Cvi  denotes  the  specific  heat  capacity  at  constant 
volume.  The  interface  exchange  terms  rk  satisfy  Y. ]krk  =  0  and  yek 
represents  molecular  and  turbulent  diffusion.  The  rate  of  energy 
deposition  Qj  results  from  fission  and  is  discussed  in  the  next 
subsection. 

The  gas  phase’s  thermal  energy  ( eg )  accounts  for  the  radiolytic 
gas  bubble  temperature.  Its  equation  is  formed  using  the  thermal 
energy  equation  for  a  bubble  as  given  in  Soo  (1990).  By  summing 
over  the  bubbles  within  a  fluid  volume  the  resulting  equation  of 
the  gas  phase  can  be  expressed  as: 

agpg^(eg)  +  agpgV  ■  ( vgeg )  =  V  •  ugyegS7Tg  +  rg, 

where  a  turbulent  diffusion  term  involving  y|  has  been  introduced 
into  this  equation. 
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3.3.  FETCH  interface  modules 

The  FETCH  interface  modules  that  link  the  fluid  and  radiation 
fields  are  described  here. 

3.3. 1.  Delayed  neutron  precursors 

The  transport  of  the  delayed  neutron  precursors  are  modelled 
by  FLUIDITY  within  the  liquid  phase  1  of  the  solution.  The  precursor 
concentration  for  delayed  group  j  is  denoted  by  a, D/r,  E,  t)  and  is 
represented  by  the  equation: 

=  _  v  ■  v,a,Dj  -  XjO (,Dj  +  Pj  J  d£lj°°  dEvaf(E)  0(fi,  E) , 

in  which  v\  is  the  velocity  of  the  liquid.  The  final  term  in  this  equa¬ 
tion  relates  to  the  source  of  precursors  and  this  is  linked  with  the 
radiation  field.  It  is  created  from  a  fraction  of  the  fission  events  that 
occur  within  the  fluid  element  where,  as  already  described,  this 
fraction  is  denoted  by  fy.  The  sink  term  comes  from  the  decay  of 
the  precursor  (with  half  life  /.j),  which,  in  turn,  will  act  as  a  source 
in  the  neutron  field.  Note  that  the  neutrons  created  from  the  decay 
of  a  precursor  will  have  their  own  unique  energy  spectrum. 

3.3.2.  Heat  transfer 

The  energy  Qj  deposited  in  the  fluid  enters  only  the  liquid  phase 
of  Eq.  (19)  and  this  arises  from  the  fission  events  within  the  radia¬ 
tion  field.  The  source  term  is  therefore  defined  as  the  reactor’s 
power  given  in  Eq.  (14): 

Qj(r,  t)=  J  dE  waf(r,E,  t)  4{r,  E.  t), 

which  is  the  total  number  of  fissions  occurring  in  the  fluid  element 
multiplied  by  the  recoverable  heat  energy  per  fission  w,  which  in 
this  case  is  taken  as  200  Mev. 

3.3.3.  Macroscopic  cross-section  nuclear  data 

The  macroscopic  cross-section  material  data  used  in  the  neu- 
tronic  calculation  is  dependent  on  both  the  solution  temperature 
and  its  gas  saturation.  The  material's  temperature  dependency  is 
resolved  by  the  pre-generation  of  tables  of  cross-sections  in 
WIMS9  (Newton  and  Hutton,  2002),  which  range  in  temperature 
from  20  °C  to  100  °C  at  10  °C  intervals  (where  it  is  assumed  the 
solution  is  pure  liquid,  i.e.  oq  =  1 ).  The  cross  sections  are  then  calcu¬ 
lated  individually  for  each  element  of  the  neutronic  mesh.  Each 
cross  section  is  found  by  linearly  interpolating  across  the  cross- 
section  tables  at  the  element’s  solution  temperature.  In  addition 
to  this,  the  density  of  the  fluid  is  then  taken  into  account  by  mul¬ 
tiplying  the  macroscopic  cross-sections  by  liquid  volume  fraction. 
This  takes  into  account  the  voidage  introduced  by  the  gas  bubble 
production. 

3.3.4.  Radiolytic  gas  production 

Whilst  the  continuity  of  the  gas  and  liquid  phases  of  the  solu¬ 
tion  are  represented  in  Eq.  (17)  an  additional  equation  for  the  dis¬ 
solved  radiolytic  gases  is  given  by: 

+  V  ■  qkockCk  =  V  ■  +  SJ  +  Fk,  (20) 

for  which  the  equation  is  solved  for  the  scalar  concentration  Q  of 
dissolved  radiolytic  gas  within  the  liquid  (k  =  /).  In  this  the  interface 
exchange  obeys  JfkFck  =  0  (described  below)  and  ffk  represents 
molecular  and  turbulent  diffusion.  The  source  for  the  production 
of  radiolytic  gas  (from  fission)  that  is  dissolved  in  fissile  liquid  is  gi¬ 
ven  by: 

$  =  J  dE  A<r'f(r,E,t)  </>{r,E,t), 

where  =  10  17  represents  grams  of  gas  liberated  per  fission 
(which  was  found  from  experiment  Barbry  and  Fouillaud  (1996)). 


The  term  oj  «  adjusts  the  fission  cross-section  so  that  it  corre¬ 
sponds  to  the  volume  fraction  of  the  liquid  phase  of  the  solution. 
The  terms  G;,  Cg  introduced  in  Eq.  (18)  can  now  be  defined  as 
Gi  =  Fj  and  Gg  =  Fg,  and  thus  the  sources  of  mass  in  the  continuity 
Eq.  (17)  is  given  by: 

S?  —  (Gi  ~  Gg),Sg  =  —(Ci  —  Cg). 

The  expression  for  the  interface  exchange  Fj  and  Fg  is  described 
in  more  detail  in  Pain  et  al.  (2001a),  which  correlates  the  rate  of 
dissolved  gas  forming  into  bubbles  in  terms  of  Henry’s  law.  How¬ 
ever,  in  fissile  solutions  a  number  of  processes  contribute  to  the 
formation  and  destruction  of  bubbles  which  include;  turbulence, 
coalescence  and  bursting  when  pearced  by  fission  fragments.  With 
so  many  unknowns,  the  modelling  of  bubble  formation  and  con¬ 
centration  is  made  to  match  experimental  observations.  In  this 
article  the  exchanges  between  dissolved  gas  and  bubbles  is  given 
as: 

Ff  =  X(C,  C0),f*  =  -Ff,  (21) 

where  X  denotes  a  bulk  diffusion  parameter  (Mather  et  al.,  1996, 
Mather  and  Barbry,  1991).  By  correlating  with  experimental  data 
and  following  closely  to  (Mather  et  al.,  1994, 1996,  Mather  and  Bar¬ 
bry,  1991),  the  following  expression  for  the  diffusion  parameter  is 
given  as: 

X  =  h(C,  -  Co)a,X(C(). 

The  value  X  is  correlated  to  match  the  pressure  traces  of  exper¬ 
imental  results  using  X(Ci)  =  N(C/  -  C0)n  in  which  C0  is  the  gas  dif¬ 
fusion  threshold  at  ambient  temperature  and  zero  solution 
pressure.  That  is,  the  pressure  in  the  bubbles  is  from  surface  ten¬ 
sion  only  when  the  bubble  has  radius  r0  =  2  x  10  2  cm.  h  denotes 
the  heavyside  function  and  its  use  prevents  gas  bubbles  from  dis¬ 
solving  back  into  the  liquid.  The  values  of  N  and  n  have  been  esti¬ 
mated  from  the  Silene  experimental  data  which  found  values 
N  =  4.71  x  1016  and  n  =  3.  It  should  be  noted  that  the  Silene  exper¬ 
iment  used  a  uranyl  nitrate  solution  and  so  is  similar  to  that  of  the 
SUPO’s  fuel. 

4.  Cooling  coil  modelling 

This  section  presents  the  formulation  for  resolving  the  heat  ex¬ 
changes  between  the  fissile  solution  and  cooling  coil  system  whilst 
the  cooling  fluid  is  pumped  through  it.  With  the  SUPO  model  rep¬ 
resented  in  RZ  geometry  it  is  not  appropriate  to  model  the  cooling 
water  system  explicitly.  Instead,  a  ID  sub-model  approach  is  used 
to  resolve  the  cooling  water’s  temperature  and  simulate  the  heat 
flow  as  the  cooling  water  advects  along  the  length  of  the  pipe.  Heat 
is  then  exchanged  between  the  ID  sub-model  and  the  SUPO  model 
by  calculating  the  appropriate  heat  transfers  on  the  regions  where 
the  pipe  surfaces  come  into  contact  with  the  fissile  solution.  It  was 
shown  in  the  previous  section  that  the  meshing  of  the  SUPO  model 
had  explicitly  resolved  the  pipe  regions  which  cut  the  RZ  plane. 
The  heat  that  exchanges  between  the  fluids  and  coil  will  therefore 
only  be  on  the  elements  that  are  adjacent  to  the  elements  occupied 
by  the  pipe.  This  section  describes  how  the  appropriate  heat  trans¬ 
fers  takes  place  together  with  the  sub-model  for  advecting  the  heat 
through  the  pipe  system. 

4.1.  Thermal  energy  equations  for  the  fissile  solution  and  pipes 

The  multi-phase  equation  for  the  temperature  of  the  fissile  li¬ 
quid  phase  T,  is  given  in  Eq.  (19).  To  develop  the  pipe  submodel, 
it  is  only  necessary  to  consider  the  advection  component  of  the 
equation,  together  with  the  heat  source  and  the  heat  exchange 
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terms  with  the  pipe.  The  temperature  equation  for  the  liquid  phase 
is  therefore  given  by: 

^  h  (J |T|  —  <7 {T coil  =  Syi,  (22) 


+  (23) 

with  Si  =  STi  +  <riTcoU  and  °t  =  u  •  V,  where  u  is  the  velocity  of  the 
liquid.  The  terms  Cpi,pi  and  ai  denote  the  fissile  liquid’s  specific  heat 
capacity,  density  and  liquid  fraction  (the  fraction  is  divided  between 
the  contributions  of  the  liquid  and  gas  phases  occupying  the  fluid 
element).  The  terms  <r,  denotes  a  sink  term  of  temperature  which 
accounts  for  the  heat  loss  from  the  fissile  solution  to  the  cooling  coil 
system,  the  coil’s  temperature  is  denoted  by  Tcoii.  Sn  denotes  the 
external  sources  of  heat  supplied  to  the  fluid.  This  comprises  of 
the  energy  generated  through  fission;  however,  in  Eq.  (23)  the 
source  term  is  also  wrapped  within  the  fluid  heat  sources/sinks 
relating  to  the  heat  exchanges  with  the  coils. 

The  method  of  discretisation  of  the  spatial  dependence  of  the 
temperature  equation  is  through  a  control  volume  method  which 
uses  suitable  flux  limiters  to  ensure  good  solution  accuracy.  The 
control  volumes  for  the  fissile  solution  temperature  are  defined 
as  the  same  finite  elements  used  to  discretise  the  domain,  and  their 


respective  function  are  denoted  by  Me,  e  e  (1,2 . /V,.lc.j,  where  Nck. 

denotes  the  number  of  control  volumes  or  elements.  These  control 
volume  functions  are  of  the  value  1  within  their  control  volume 
and  0  outside,  thus  the  resulting  temperature  profile  will  be 
piece-wise  constant  over  the  control  volumes.  This  approximation 
is  given  by: 

T,  =  J2MJ,e.  (24) 

The  heat  equation  that  is  solved  for  each  pipe  in  the  model  is  gi¬ 
ven  by: 

(25) 

or  alternatively, 

^  +  Cr.p,^  +  «.T,.S„  (26) 

where  Sw  =  awTt  (the  over  line  indicates  that  the  heat  source  comes 
from  the  main  FETCH  model  and  this  will  be  averaged  in  some  sense 
as  described  below).  These  equations  resolve  the  temperature  of  the 
coil’s  water  Tw  and  steel  Tst  components.  The  terms  Cp  and  p  denote 
the  steel  and  water  specific  heat  capacity  and  density,  respectively, 
and  VR  denotes  the  volume  ratio  of  the  pipe’s  steel  and  cooling 
water.  Similar  to  the  previous  equation,  the  term  nw  denotes  the 
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Fig.  7.  The  real  time  series  simulation  of:  (a  and  b)  the  maximum  gas  and  liquid  velocities;  (c  and  d)  the  maximum  gas  and  liquid  temperatures. 


Fig.  8.  The  graph  show  the  temperature  distribution  of  the  cooling  water  along  the 
three  pipes  whilst  SUPO  is  operating  at  a  quasi-steady-state  25  kW  power.  The  lines 
coloured  black,  red  and  green  represent  the  distributions  of  the  first,  second  and 
third  cooling  coil  sets  respectively.  (For  interpretation  of  the  references  to  colour  in 
this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 

sink  (or  in  this  case  source  with  T;  >  Tw )  terms  that  represents  the 
heat  exchanging  from  the  fissile  solution  to  the  pipe.  These  source 


and  sink  terms  for  exchanging  heat  are  described  in  more  detail 
below. 

In  resolving  the  time  dependence  of  Eq.  (26)  sub-cycling  is  used 
with  time  step  sizes  A tpipe  =  min  j  At,  0.1  A  high-resolution 

control  volume  finite  element  method  is  applied  to  discretise  the 
heat  equation  in  space.  In  this  approach  the  spatial  domain  is 
divided  into  a  finite  number  of  non-overlapping  elements  over 
which  the  heat  varies  with  a  piece-wise  constant  variation.  That 
is,  the  temperature  of  the  pipe  water  is  represented  as: 

Tw  =  J>Twe,  (27) 

where  CL  are  the  control  volume  basis  functions  (which  have  the 
value  1  inside  their  element  and  0  elsewhere)  and  Tw  e  their  corre¬ 
sponding  expansion  coefficients. 

4.2.  Calculating  the  heat  transfer  coefficient 

The  transfer  of  heat  from  the  fissile  solution  to  the  pipe  is  cor¬ 
related  from  the  Nusselt  number: 

Nu  =  th  +  a2(RePipePrfluid)^,  (28) 

for  which  the  values  used  are  a^  =  5  and  a2  =  0.025.  Using  SGI  units 
(cm  and  g)  the  values  of  the  Prandtl  number  is  given  as  Ptyuid  = 
2.842,  the  kinematic  viscosity  of  the  fissile  solution  is  given  by 
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Fig.  9.  Time-averaged  of:  (a)  gas  volume  fraction  and  (b)  fissile  solution  temperature  (in  °C). 


Fig.  10.  Top  and  bottom  sections  of  the  SUPO  reactor  showing  the  time  averaged  liquid  (a  and  b)  and  gas  (c  and  d)  phase  velocity. 


Mfiuid  =  0.004322  and  the  thermal  conductivity  kfluid  =  0.00656.  The 
pipe  Reynolds  number  is  Repipe  =  Taking  into  consideration 

the  scaling  for  the  volume  fraction ^of  the  fissile  solution  (a;),  the 
heat  transfer  coefficient  is  given  by: 


(29) 


In  practice  the  heat  transfer  coefficient  will  be  calculated  for 
each  element  e  on  the  fissile  solution  mesh.  This  is  given  as: 


h  _  a,e(Nu)e  kjm 

<W 


(30) 


where  (Nu)e  denotes  the  average  value  of  Nu  inside  element  e.  These 
values  are  used  in  the  parameterisation  of  heat  transfer  at  the  unre¬ 
solved  scales  for  which  the  rate  of  exchange  of  heat  is  given  by: 


y  =  h(Tw-T(),  (31) 

where  n  is  the  unit  vector  which  is  normal  to  the  pipe  and  points 
inwards.  Multiplying  Eq.  (26)  by  the  control  volume  basis  function 
Qi  and  integrating  over  the  pipe  one  obtains: 

Jv  Q.,(VRCPstpJTd»  +  CPwPwDu^)dV+J  Qi(hT„-(hTi))dr=0,  (32) 

where  Vpipe  is  the  pipe  volume  and  rpipe  the  pipe  surface.  By  assum¬ 
ing  the  water  and  steel  to  have  the  same  temperature,  Tst  =  Tw,  Eq. 
(32)  is  simplified  to: 

^  +  G„a.  (^r)i  +  ^  (hr., -  (B»,)  -  0, 

(33) 
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Fig.  11.  (a)  Time-averaged  shortest-lived  delayed  neutron  precursor  concentration  and  (b)  dissolved  radiolytic  gas  concentration  (g  cm  3). 


Fig.  12.  These  pictures  show  snap  shots  of  SUPO  at  time  instance  2264  s  into  the  transient  simulation,  (a)  Gas  volume  fraction,  (b)  temperature,  (c)  longest  delayed  neutron 
precursor,  (d)  shortest  lived  delayed  neutron  precursor  (e)  concentration  of  dissolved  gas. 


where  Cpjpe  outer,  Awater  and  Aouter  denote  the  pipe’s  outer  circumfer¬ 
ence,  cross  sectional  area  occupied  by  the  water  and  total  cross 
section  area.  The  top  hat  marks  the  heat  transfer  variables  that 
are  mapped  from  the  FETCH  model  and  these  will  be  described  in 
the  next  subsection.  VR  -  By  comparing  Eq.  (33)  with 

Eq.  (26)  the  absorption  <jw  and  the  source  Sw  values  on  pipe  at 
control  volume  i  can  now  be  estimated  as: 


Swi  =  C)je  outer  (ftT)),-  (35) 

4.3.  Heat  exchanges  between  the  pipe  sub-model  and  the  RZ  FETCH 
SUPO  model 

In  order  to  transfer  information  between  a  cooling  coil  and  the 
RZ  FETCH  model,  a  set  of  basis  functions  Gk(s),  fce{l,2,...,JVcoa}, 
are  defined  where  JVcoil  denotes  the  number  of  intersections,  or 
turns,  that  the  coil  makes  with  the  RZ  plane.  The  variable  s  denotes 
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the  arc  length  along  the  pipe  and  so  varies  from  sin  =  0  at  the  pipe’s 
inlet  to  Sout  at  its  outlet.  Each  function  Gfc(s)  is  associated  with  the 
kth  intersection  point  and  centres  over  a  node  positioned  at  sk 
which  is  the  distance  along  the  pipe  where  the  kth  intersection  oc¬ 
curs.  These  functions  are  similar  to  standard  ID  finite  elements  in 
that  they  have  the  value  1  over  their  respective  node  and  fall  away 
linearly  to  the  value  zero  on  their  neighbouring  nodes.  Beyond 
their  neighbouring  nodes,  the  functions  have  zero  value. 

In  transferring  heat  between  the  models  the  following  lumped 
variables  are  defined: 

Lk  =  f  Cuds,  (36) 

which  essentially  provides  a  measure  of  the  length  of  pipe  about 
intersection  k,  and 


T  fcoilGkTwds 

Wk  fcouCkds  : 


(37) 


which  is  the  coil’s  water  temperature,  weighted  with  the  basis  func¬ 
tions  Gfc,  and  integrated  across  the  length  of  the  coil.  This  provides 
an  average  temperature  of  the  pipe’s  water  about  intersection  k. 
In  addition,  the  heat  transfer  coefficient  expressed  at  the  centre  of 
the  coil’s  control  volume  Mf  (this  is  at  point  point  s  say)  is  given  by: 


h(s)  =  5>(s)fik. 


(38) 


This  expression  is  used  to  define  the  absorption  term  of  Eq.  (34). 
Its  value  is  dependent  on  the  average  heat  transfer  coefficient  from 
the  fluid  mesh,  and  this  is  calculated  through  a  volume  average  of 
the  heat  transfer  coefficient  about  the  region  of  intersection  k  (in 
the  RZ  model).  That  is,  at  each  intersection  k  the  heat  transfer  coef¬ 
ficient  is  calculated  as: 


c  Ee^IyMehedV 

k  EeeSkfvM'dV  ' 


(39) 


where  the  summation  runs  over  the  elements  of  the  fluids  mesh 
that  contain  fissile  solution  and  that  are  adjacent  to  the  pipe’s  kth 
intersection  (this  set  is  denoted  by  £k).  A  similar  expression  is  used 
to  generate  an  average  heat  exchange  term  centred  about  the  kth 
coil  intersection: 


(hTl)k  = 


Eeeek  JvMeheTledV 
Eee  ekJvMedV  ’ 


(40) 


which  can  then  be  used  to  represent  the  heat  exchange  that  takes 
place  on  the  ID  pipe  sub-model  over  control  volume  i: 


(hTi)i^'EG*(sd(hTl)k. 


(41) 


This  term  provides  the  required  expression  for  the  source  term 
in  Eq.  (35)  to  be  calculated. 

It  now  remains  to  calculate  the  absorption  term  of  the  heat 
equation  in  the  RZ  SUPO  model  that  accounts  for  the  heat  lost  to 
the  cooling  coils.  The  approach  is  to  take  the  temperature  Eq. 
(23)  for  the  fissile  liquid  and  to  multiply  it  by  a  control  volume  ba¬ 
sis  function  Me,  and  integrate  this  over  the  domain  V.  This  gives  the 
following  expression: 

JvMeCPlaiPl(^Pj  dV+ Mseh(T,-Twm)dr 
-  J  MeSndV,  (42) 

where  (.)e  expresses  the  value  of  the  variable  at  the  centre  of  ele¬ 
ment  e  and  Twfc(e)  is  the  average  pipe  temperature  about  an  intersec¬ 
tion  K(e)  -  see  Eq.  (37).  The  term  Mse  is  given  by: 


J 


JyMedV r 

We  p 


outer  Lm, 


(43) 


where  We  =  ]Trar£.  f  Mm  dV.  The  index  set  K(e)  =  k  denotes  the 
pipe  intersection  k  that  is  adjacent  to  the  fluid  element  e.  The  sum¬ 
mation  in  We  is  therefore  a  volume  sum  of  all  the  fluid  elements 
that  surround  the  intersection  k.  This  is  used  to  divide  the  volume 
of  fluid  element  e  ( JvMedV )  to  give  its  relative  volume  fraction  of 
all  the  elements  that  surround  the  intersection.  This  fraction  is  then 
multiplied  by  the  total  surface  area  of  the  pipe  that  surrounds  the 
intersection  k  ( CPipe  outer  ix(e))  so  that  Mse  gives  an  expression  of 
the  surface  area  of  the  pipe  in  contact  with  element  e. 

Eq.  (42)  is  now  divided  by  the  volume  of  the  element  e  to  give: 


(DUlT,\  he(Tie  -  Twm) 


fvMedV 


^  J  Msedr  =  sTle. 


(44) 


By  comparing  the  terms  in  Eqs.  (44)  and  (23),  the  temperature 
absorption  terms  for  the  fissile  liquid  phase  can  be  expressed  as: 


Ole 


he  C„iPe  outer  JcoilGmdS 

We 


(45) 


in  which  he  is  the  heat  transfer  coefficient  (Eq.  (29))  in  element  e. 
The  temperature  source  for  the  liquid  (fissile  solution)  phase  for 
an  element  e  adjacent  to  a  pipe  intersection  is  given  by: 

Sle  =  aleTwm+STle.  (46) 


5.  Results  from  the  default  configuration  25  kW  power  SUPO 
reactor  with  initial  solution  temperature  60  °C 

The  FETCH  simulations  presented  here  are  based  on  SUPO  oper¬ 
ating  at  the  intended  25  kW  of  power.  The  aim  is  not  to  follow  a 
power  rise  to  operating  conditions  in  full  detail,  but  to  begin  with 
initial  conditions  which  should  lead  to  quasi-equilibrium  reason¬ 
ably  quickly.  The  simulations  begin  with  the  reactor  having  a  solu¬ 
tion  temperature  of  60  °C  and  a  5  °C  cooling  water  at  the  pipe’s 
inlet.  The  simulations  are  initialised  with  the  control  rod  position 
chosen  so  that  the  reactor  is  in  an  exactly  critical  state.  Each  sim¬ 
ulation  is  then  allowed  to  evolve  where  by  the  control  rod  adjusts 
its  position  over  time  in  response  to  the  reactor’s  power  so  that  it 
tries  to  maintain  25  kW  of  power.  The  heat  losses  from  the  fluid, 
other  than  through  radiolytic  gas  release  (which  are  small),  are  as¬ 
sumed  come  only  from  the  cooling  water  system.  There  will  there¬ 
fore  be  a  balance  between  the  heat  energy  extracted  by  the  cooling 
water  and  the  fission  power  once  a  reactor  reaches  a  quasi-steady- 
state.  It  should  be  noted  that  due  to  the  nature  of  multiphase  flow 
and  the  various  feedback  mechanism,  this  system  is  not  expected 
to  reach  a  true  steady  from,  but  instead  oscillate  about  some  aver¬ 
age  state. 

In  total,  four  simulations  have  been  performed  which  vary  in 
the  speed  in  which  their  control  rod  can  move  in  order  to  keep 
the  reactor  operating  at  its  intended  power.  The  first  simulation 
used  a  control  rod  that  moved  with  velocity  0.01  cm/s,  which 
should  only  be  enough  to  respond  to  long  term  power  trends. 
The  second  simulation  used  a  control  rod  velocity  of  1  cm/s. 
Although  this  is  considered  to  be  quite  excessive  in  a  reality,  it 
has  been  used  here  to  demonstrate  its  ability  to  keep  up  with  the 
short  term  power  trends  including  those  associated  with  transient 
bubble  formation.  In  order  to  eliminate  the  effect  of  the  control  rod 
movement  a  third  simulation  was  performed  with  the  control  rod 
set  in  a  static  position.  In  this  case  the  rod’s  height  was  chosen  so 
that  the  reactor  achieved  on  average  approximately  25  kW.  Finally 
the  fourth  simulation  presents  a  mesh-convergence  test.  This  sim¬ 
ulation  used  a  rod  velocity  0.01  cm/s,  but  the  spatial  discretisation 
was  increased  in  resolution  to  that  used  in  the  first  simulation.  In 
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all  calculations  the  SUPO  reactor  was  simulated  for  40  min  of  real 
time  from  its  initial  condition  stated  above. 


5.1.  The  initial  SUPO  simulation 

This  section  presents  the  first  simulation  of  the  SUPO  reactor 
that  was  allowed  to  evolve  over  time  whilst  being  regulated  by 
the  control  rod  with  velocity  0.01  cm/s.  The  graphs  presented  in 
Fig.  6  show  the  real  time  series  of  SUPO’s  maximum  power  and 
pressure,  mean  outflow  of  cooling  water  temperature  and  heat 
transfer  coefficient.  The  power  shown  in  Fig.  6a  is  remarkably  stea¬ 
dy  and  changes  by  about  5%  after  the  initial  transient.  This  result 
highlights  the  fact  that  the  SUPO  simulation  does  not  reach  a  stea¬ 
dy  state.  Instead  it  has  formed  into  a  quasi-steady  state  since  the 
power  trace  exhibits  a  noise  due  to  the  various  feedback  effects  - 
eg.  gas  bubble  production  and  its  eventual  escape  through  the 
solution  surface.  Fig.  6b  shows  the  maximum  pressure  which  again 
shows  a  steady  profile  once  the  reactor  has  reached  its  quasi-stea- 
dy-state.  The  remaining  graphs  of  Fig.  6  show  the  mean  cooling 
water  outlet  temperature  to  be  33  °C  and  the  mean  heat  transfer 
coefficient  to  be  0.14  W/cm/K.  Again  these  two  profiles  are  subject 
to  noise  and  this  would  be  a  result  of  the  fluctuations  observed  in 
the  reactor’s  power.  The  maximum  gas  and  fissile  solution  veloci¬ 
ties  and  temperatures  are  shown  in  Fig.  7.  This  shows  the  maxi¬ 
mum  fissile  solution  temperature  to  only  vary  by  a  few  degrees 
(Fig.  7d).  The  distance  along  each  of  the  three  pipes  versus  cooling 
water  temperature  is  shown  in  Fig.  8.  This  shows  a  substantial  var¬ 
iation  in  the  cooling  water  temperature  (^25-30  °C)  within  the 
pipes  for  which  the  third  (or  outer)  pipe  is  approximately  10% 
higher  than  the  two  inner  coils.  From  this  simulation  it  was  found 


that  once  the  reactor  reached  its  quasi-steady-state,  the  mean  fis¬ 
sile  solution  temperature  was  72  °C. 

It  should  be  noted  that  the  heat  removal  from  the  vessel,  as¬ 
sessed  on  the  basis  of  the  modelled  temperature  rise  across  the 
coolant  tubes,  exceeds  by  about  10%  that  derived  from  the  data 
in  King  (1956).  However  these  referenced  estimates  take  account 
of  two  additional  heat  loss  mechanisms  -  heat  loss  to  the  reactors 
recombiner  and  to  vessel’s  the  graphite  reflector.  Together  these 
amount  to  8%  of  the  total  heat  and  this,  in  turn,  will  reduce  the  de¬ 
mand  to  transfer  heat  to  the  cooling  coils.  As  the  current  FETCH 
model  assumes  all  fission  heat  appears  as  a  temperature  rise  in 
the  solution,  and  that  heat  removal  is  only  though  the  cooling  coils, 
an  increase  in  cooling  coil  temperature  would  be  expected.  The 
simulated  results  are  therefore  encouraging  and  highlight  the  addi¬ 
tional  heat  loss  to  be  an  area  for  more  detailed  analysis.  This  might 
also  explain  the  temperature  difference  between  the  predicted 
70  °C  compared  with  the  expected  60  °C. 

Figs.  9-1 1  present  the  ’time  averaged’  profiles  of  reactor’s  field 
variables  associated  with  both  the  neutronic  and  fluid  dynamics. 
By  time  average  it  is  meant  that  the  average  of  each  field  variable 
is  taken  over  the  final  30  min  of  the  simulation  -  the  time  period 
where  the  reactor  was  operating  in  its  quasi-steady-state.  The 
time-averaged  gas  volume  fraction  and  temperatures  are  pre¬ 
sented  (Fig.  9)  and  show  a  gradual  increase  in  gas  volume  fractions 
with  height.  There  are  also  low  and  high  gas  volume  fractions  on 
above  and  below  the  cooling  coils,  respectively,  and  this  shows 
the  gas  becoming  trapped  and  accumulating  under  them.  The  tem¬ 
perature  of  the  solution  outside  the  cooling  coils  varies  only  mar¬ 
ginally  by  a  few  degrees. 

The  time  averaged  fissile  liquid  and  gas  velocities  are  shown  in 
Fig.  10.  It  can  be  seen  clearly  that  there  are  two  distinct  regions  of 
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Fig.  14.  Time  series  of  several  fields:  (a)  power  (in  W),  (b)  average  solution  temperature  (in  °C)  (c)  mean  outlet  cooling  water  temperature  (in  °C)  and  (d)  mean  heat  transfer 
coefficient  in  the  cooling  water  -  fissile  solution  system  (in  W  crrr1  K  ').  Black  is  the  fine  mesh  standard  case.  Red  the  simulation  result  with  a  fixed  control  rod.  Green  is  with 
slow  control  rod  movement  and  blue  with  fast  control  rod  movement.  (For  interpretation  of  the  references  to  colour  in  this  figure  legend,  the  reader  is  referred  to  the  web 
version  of  this  article.) 


flow.  In  the  top  region  of  the  reactor  the  liquid  moves  down  in  the 
central  region  and  up  at  the  sides.  In  the  lower  region  the  liquid 
rises  in  the  vessel’s  centre  (due  to  buoyancy  forces  generated  by 
radiolytic  gas  bubbles)  and  falls  near  the  walls.  This  results  in 
two  distinct  regions  of  isolated  flow.  The  shortest-lived  delayed 
neutron  precursor  concentration  is  shown  in  Fig.  11a,  which  illus¬ 
trates  the  power  distribution  of  SUPO.  The  power  is  relatively  low 
in  the  upper  region  and  this  is  due  to  the  presence  of  the  control 
rod.  The  time-averaged  concentration  of  dissolved  radiolytic  gases 
(see  Fig.  lib)  is  fairly  uniform  through  the  solution  and  is  approx¬ 
imately  equal  to  the  critical  concentration  (given  by  Henry’s  law) 
of  radiolytic  gases  dissolved  in  the  solution  at  which  bubbles  start 
to  form. 

A  snap  shot  of  SUPO  taken  2264  s  into  the  simulation  is  pre¬ 
sented  in  Fig.  12.  At  this  time  instance  the  quasi-steady-state  con¬ 
ditions  has  been  attained  and  Fig.  12a  shows  the  gas  volume 
fraction  to  be  highly  non-uniform.  However  the  temperature  is  still 
fairly  uniform  as  shown  in  Fig.  12b.  Non-uniformity  is  also  seen  in 
the  longest-lived  delayed  neutron  concentration  field  that  is 
shown  in  Fig.  12c.  However,  the  shortest-lived  delayed  concentra¬ 
tion  (or  power  distribution)  in  Fig.  12d  looks  similar  to  the  time- 
averaged  fields  shown  in  Fig.  11a.  Also,  the  concentration  of 
dissolved  gas  is  still  fairly  uniform  and  near  the  critical  concentra¬ 


tion,  this  is  shown  in  Fig.  12e.  The  instantaneous  velocities  shown 
in  Fig.  13  still  show  the  two  circulation  regions  mentioned  in  the 
time  average  distributions. 

5.2.  Analysis  of  the  transient  simulation  with  variation  in  the  spatial 
discretisation  and  control  rod  speed 

The  numerical  simulations  of  the  four  test-cases  with  varying 
control  rod  speed  and  spatial  meshes  are  compared.  Fig.  14  shows 
the  time-series  of  several  of  the  reactor’s  fields  including:  power, 
mean  solution  temperature,  mean  outflow  cooling  water  tempera¬ 
ture  and  mean  heat  transfer  coefficient.  The  power  versus  time 
plot,  Fig.  14a,  shows  that  the  variation  in  fission  rate  is  minimal 
with  the  fastest  control  rod  and  varies  by  approximately  15%  for 
the  static  control  rod.  Again  this  shows  that  the  reactor  enters  a 
quasi-steady-state  for  all  simulations.  They  all  oscillate  about  some 
mean  state  with  the  noise,  due  to  the  feed  back  effects,  being 
dampened  with  a  faster  responding  control  rod.  Average  solution 
temperature  for  all  cases  (Fig.  14b),  has  a  small  variation  of 
0.2  °C  and  the  mean  outlet  cooling  water  temperature  versus  time 
(Fig.  14c)  varies  by  1  °C  for  all  cases.  The  mean  heat  transfer  coef¬ 
ficient  (averaged  along  all  the  pipes)  (Fig.  14d),  changes  by  approx¬ 
imately  5%  over  the  course  of  the  simulation. 
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Fig.  15.  Time  series  of  control  rod  tip  position  verses  time  (in  cm  from  the  vessel's 
base).  Black  indicates  the  fast  moving  rod  on  the  fine  mesh,  red  indicates  the  fixed 
rod  on  the  standard  mesh,  green  the  slow  moving  rod  on  the  standard  mesh  and 
blue  the  fast  rod  on  the  standard  mesh.  (For  interpretation  of  the  references  to 
colour  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


The  time-series  of  control  rod  tip’s  position  (from  the  vessel’s 
base)  versus  time  is  shown  in  Fig.  15.  These  show  the  control  rods 
continually  moving  in  and  out  of  the  vessel  which  will  be  in 
response  to  the  noisy  power  profile.  The  simulation  with  rapid 
1  cm/s  rod  movement  shows  a  large  variations  in  the  tip  position 
while  the  slow  movement  case  (0.01  cm/s)  shows  little  variation 
as  the  response  is  effectively  time  averaged.  Finally,  using  the  finer 
spatial  mesh  serves  to  decrease  the  simulated  reactivity  slightly, 
and  this  results  in  a  higher  control  rod  tip  position  but  with  a  sim¬ 
ilar  variation  about  its  mean. 

6.  Conclusion 

The  Modelling  of  the  time  dependent,  coupled  neutronic  and 
fluid  dynamics  of  the  SUPO  reactor  has  been  conducted  using  the 
FETCH  model.  The  simulation  has  modelled  the  reactor  in  RZ 
geometry  which  resolves  the  fuel  solution  in  terms  of  multi-phase 
flow  (consisting  of  both  liquid  and  gas)  in  conjunction  with  the 
neutronics,  together  will  their  feedback  effects.  The  models  are 
run  with  a  semi-explicit  representation  of  the  cooling  coils  and  a 
’sub  model’  method,  which  has  been  developed  here,  is  then  used 
to  transfer  heat  from  the  bulk  fluid  into  the  coil’s  water.  An  adap¬ 
tive  control  rod  movement  algorithm  with  varying  velocity  is  also 
used  in  the  simulation.  Using  this  algorithm  the  rod’s  insertion 
length  is  constantly  adjusted  to  respond  to  changes  in  reactivity 
with  the  aim  to  keep  the  reactor  operating  at  25  kW.  Overall, 
FETCH  is  an  excellent  tool  to  analyse  dynamics  of  SUPO.  This  inves¬ 
tigation  has  not  only  shed  some  light  on  its  dynamics  that  were  not 
previously  recorded  during  its  operation  (such  as  flow  patterns  and 
long  term  power  trends),  but  also  the  simulation  agreed  well  with 
data  that  was  recored,  such  as  its  heat  transfer.  The  conclusions  of 
this  study  are: 

(a)  When  modelling  the  25  kW  condition  data  set  using  the 
adaptive  control  rod  algorithm  to  achieve  quasi-steady 
power  conditions,  it  has  been  assumed  in  the  FETCH  model 
that  all  fission  heat  arising  in  the  solution  is  removed  by 
the  three  cooling  coils.  The  predictions  of  the  modelled  cool¬ 
ant  outlet  temperatures  and  those  measured  in  SUPO  itself 
are  encouragingly  close  when  the  estimated  6%  heat  loss 
by  conduction  to  the  surrounding  graphite  reflector  is  taken 
into  consideration.  This  loss  will  be  accounted  for  in  future 


modelling  so  that  any  residual  differences  in  solution  tem¬ 
perature  may  be  ascribed  to  errors  in  solution-to-coils  heat 
transfer  algorithms  (taken  together  with  any  effects  of  mov¬ 
ing  from  2D  to  3D  FETCH  modelling).  We  consider  that  both 
heat  losses  through  conduction  and  losses  due  to  escape  of 
the  ’mobile’  part  of  the  fission  energy  release  residing  in  fast 
neutrons  and  gamma  rays  (bearing  in  mind  the  very  small 
SUPO  core)  should  be  assessed. 

(b)  Using  two  different  averaging  times,  modelled  high  resolu¬ 
tion  traces  of  control  rod  movement  in  maintaining  a 
quasi-steady  power  have  been  obtained.  The  observed  fine 
scale  time  structure  is  due  to  the  reactivity  feedback  effects 
of  bubble  formation,  movement  etc  and  also  its  effect  on  free 
surface  levels.  This  fine  scale  trace  could  be  analysed  to  yield 
further  information  on  the  stability  of  SUPO  at  its  normal 
operating  power  -  or  at  higher  powers  which  were  utilised 
earlier  in  the  reactor  life  and  when  suggestions  were  made 
regarding  localised  boiling  between  coils. 

(c)  There  are  two  basic  flow  recirculation  zones  that  form  in  the 
2D  simulations.  It  is  not  clear  if  this  would  be  the  case  with 
3D  modelling  of  SUPO.  SUPO  modelling  shows  a  relatively 
homogeneous  temperature  distribution  with  the  tempera¬ 
ture  varying  only  a  few  degrees  across  the  reactor.  There  is 
a  gradual  increase  in  time-averaged  gas  volume  fraction  as 
we  move  up  through  the  solution.  Remarkably,  this  average 
is  (with  the  exception  of  around  the  cooling  coils)  approxi¬ 
mately  uniform  in  the  horizontal. 

(d)  We  notice  in  the  modelling  a  fairly  dense  liquid  with  few 
bubbles  immediately  above  individual  cooling  coils  -  with 
a  relatively  large  liquid  volume  fraction,  in  a  time  averaged 
sense,  below  the  tube.  Below  a  tube  there  is  a  concentration 
of  trapped  gas  whereas  above  it  there  is  a  relatively  low  con¬ 
centration  of  radiolytic  gases  due  to  the  gas  bubbles  leaving 
this  area  and  not  being  replace  from  below. 
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